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1. Introduction 

Lattice models have long been a powerful tool to study quantum mechanics and 
quantum field theories. They allow for computational treatment of analytically 
intractable problems, from high-Tc superconductors [1 to colour confinement in 
QCD [gj- However, the computational cost of dealing with these problems can 
be exorbitant [3]. In the spirit of Feynman [4 , great effort has been devoted to 
simulating these models in intrinsically quantum systems, which can be more efficient 
in reproducing key characteristics of these problems. Indeed quantum simulators 
have managed to produce results in the study of relativistic quantum mechanics and 
quantum phase transitions O |6] . 

Optical lattices are one of the most promising candidates for quantum simulations, 
due to their scalability, tunability and their versatility in terms of geometry and 
dimensionality. They have already been successful in simulating quantum phase 
transitions of Bose-Hubbard and spin Hamiltonians [II [H [9] , and several proposals 
have been made to observe synthetic gauge field theories and topological order [TOl 
\TT\ [T2j [131 [H]- III this work we propose a way to study in an optical lattice the 
2+1 Dirac fields coupled to both Abelian and non-Abelian gauge fields. Our scheme 
relies on intensity modulations of the trap and of additional laser beams to give birth 
to all sorts of phenomena: from the appearance of synthetic electromagnetic fields 
and a Dirac mass to exotic fiavour-coupling perturbations. Our idea is intimately 
connected to the way mechanical deformations induce effective magnetic fields in 
graphene [15l [161 El ? but makes use of the greater tunability of optical lattices to 
create a greater family of potentials. 

This work is structured as follows. In section |2] we show how the Dirac equation 
arises from a tight-binding treatment of the honeycomb lattice. We pay special 
attention to the tools for deriving continuum limit and to the appearance of an extra 
degree of freedom, called fiavor. In section |3] we analyse the effect of perturbations on 
the lattice, interpreting the result as a coupling between the Dirac field and scalar, 
Abelian and non-Abelian external gauge fields. In section [4] we detail the experimental 
implementation of the Dirac model and of the external potentials using ultracold atoms 
in an optical lattice. After discussing a trapping scheme based on state-dependent 



optical lattices, section 4.3 shows how to implement the effective potentials from 
section [3] using intensity modulations of the laser beams. Finally, section [5] discusses 
different experimental protocols to characterize and measure the previous quantum 
simulations. 

2. From tight-binding models to the Dirac fields 

Here we develop a theoretical framework in which 2D fields arise as the continuum 
approximation to a quadratic tight-binding model on a honeycomb lattice. Later in 
section |4] we show how these Hamiltonians describe non-interacting atoms in a strong 
periodic confining potential, but in this section the focus is on the abstract model. 
More precisely, we are interested in the derivation of honeycomb lattice band structure, 
whose excitations behave as relativistic Dirac particles, and how this procedure is 
modified to obtain a coupling between the Dirac particles and emergent gauge fields. 
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Figure 1. Two triangular sub-lattices displaced one with respect to the other to 
generate a honeycomb lattice. The vectors vi and V2 are the hexagonal lattice 
generator vectors, vq = and V3 = vi -\-V2. The unit cells are shown with dashed 
rectangles, and are labelled by m. Ta and Tb are the hopping parameters related 
with jumps between positions inside sub-lattices A and B while Ji are the hopping 
parameters related with hexagonal jumps between sublattices, each one related 
with the corresponding Vi, i = 0, 1, 2. 



2.1. Tight-binding model for two coupled sub-lattices 

Consider a system of fermionic particles in a strong confining periodic potential. 
We model these particles with fermionic creation and annihilation operators which 
describe the presence or absence of a fermion in a particular lattice site. The lattice 
is bipartite, which means that it can be divided into two disjoint sets of sites A and 
B, where the nearest neighbours of any A-site are all elements of 5, and viceversa 
[cf. figure lla)]. Thanks to the implementation we may consider nearest-neighbour 
hoppings between different sublattices of a bipartite lattice, and also consider next- 
to-nearest-neighbour (or intra-lattice) hoppings. For the hexagonal spatial geometry 
shown in figure [l] our model Hamiltonian reads: 
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where the Ji are the hopping parameters related with the first neighbour hoppings; 
r^ and r^ model the tunneling effect inside sublattices A and B respectively; aj^ 
and am {bl^ and bm) are the creation and annihilation operators for an atom in the 
position m of sublattice A (5), and the v^s refer to the v vectors represented in figure 
[I]3. The first line of equation ^ represents the hexagonal hoppings, and the second 
line represents the two kinds of triangular hoppings. 



2.2. Energy bands and Dirac field 

The honeycomb Hamiltonian is recovered from equation M by making Jq 
J2 = J and Ta = Tb = 0: 
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Figure 2. Left: First Brillouin Zone (IBZ) for an hexagonal lattice. The K± 
points written in black are the non-equivalent Dirac Points. Center: three- 
dimensional representation of the energy bands of the honeycomb lattice. The 
two surfaces represent the positive and the negative energy solutions equation Cj 
meeting together in the Dirac points where E = 0. Right: zoom around the Dirac 
points where the energy bands take conical shape and the Hamiltonian of the 
system can be approximated as a Dirac Hamiltonian, being E oc \q\, equation (|6|. 



To calculate the energy bands, we do a Fourier Transform (FT) over the creation and 
annihilation operators as 

d^k 



.t 



-ik-Vj 



27r 



•at 



(3) 



where k is the momentum, Vrn is the position of the site m^ the integration limit Q. 
makes reference to the first Brillouin Zone (IBZ) and a similar expression is used for h. 
The diagonalization of the Hamiltonian in momentum space leads to the eigenvalues 
of the energy that form the energy band structure; namely 



£±(fe) = ± 



\ 



3 + 4 cos ( -/Ca; 1 COS 



■V3 
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which are represented in figure [2j There are six points in the IBZ for which £^ = 0, but 
only two of them are non-equivalent and we name them K± = (=F^, 0)/(i, where d is 
the unit cell size [18]. Assuming only small energy perturbations around £^ = 0, the 
excitations will be confined to the so called Dirac cones. To find the effective theory 
for these low-energy and long-wavelength excitations one makes an expansion around 
the Dirac points of E{k)^ introducing q = k — K and assuming that for the relevant 
states it remains "small" (i.e. \q\ x d <^ 1). Under these conditions, it is possible to 
work with the Dirac cones as if they extended over all values of g, with integrals in d^k 
over the IBZ being replaced by an integral in d^q over the whole momentum space. 
By transforming the operators back to position space via inverse FT 

^iir-) = I l-a^.^^e-^(^+'^-^^d\ (5) 

one gets an effective theory for the continuum fields {(/)a(r), (/)5(r)}, satisfying the 
Dirac-like Hamiltonian (cf. equation [14]) 

h+{q) ^a-q, h-{q) « -h+{q)* , (6) 

where a — (ox^^y) are the usual Pauli matrices and the subscript ± stands for the 
choice of cone {K± ) . 
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2. 3. Derivation of the Dirac Hamiltonian in position space 

While the previous discussion reveals the main ingredients of the relativistic fields 
which arise from our discrete Hamiltonian, two concerns appear. One is the added 
difficulty of treating spatially dependent perturbations, caused by the two consecutive 
Fourier transforms. The other concern is the little attention paid to the fact that there 
are two intrinsically different cones, which in line with previous literature we call the 
two "ffavors" of Dirac particles in the lattice. As we will show later, this ffavor plays 
a role in the simulation of non-Abelian fields. 

In order to extend the derivation of the Dirac dispersion relation to setups 
in which translational invariance is weakly broken, we rely on a continuum-field 
approach that bypasses the use of momentum space [19]. Starting with our tight- 
binding Hamiltonian (fTl) with suppressed intralattice hoppings {Ta = ^b = 0), we 
approximate the Fock operators at each lattice site as the value of a continuous field 
defined over all space, but which varies so smoothly that it is approximately constant 
over each unit cell. Moreover, since we are interested in excitations around the K± 
quasimomenta, these fields are the envelope of quasi-plane wavepackets around those 
points 

a« ^ Vd^ (^„+(r-™)e-^^+^^'" + *„_(f„)e-^^-^^'") , (7) 

with an equivalent expression for hm- Here d^r is the unit cell area and the 
exponential containing the Dirac points K± implements the desired wavepacket ansatz. 
Substituting these operators in the tight-binding model produces the following limit 
Hamiltonian 

HocY^ I e'^^^^-^-^ [^i^{r)^hv{r^Vi)e-'^^^^'\ d^r, (8) 

where the Greek indices stand for the two possible fiavour choices (± cones) and the Vi 
stand for the vectors 'iTo7'^i7'^2 depicted in figure [l] Since we have enforced the fields 
to be slowly varying, the presence of the rapidly oscillating exponential g^^C^^-^^) 
imposes the condition ji = v and we end up with two fiavour-decoupled integrals. At 
this point one may expand ^{f ^ v) c^ [1 + i7V]^(r) so the Hamiltonian reads: 

^ °^ E / *a/'(^) [l + ^i"^] *6p(^le-*^''"M2r. (9) 

Our resulting Hamiltonian is therefore diagonal both in position and in fiavour space 
but couples the internal degrees of freedom a-b via an off-diagonal operator with only 
one nontrivial term: 

C^ = 1 + e-^^'^^"^^ [1 + vW] + e'^^'"'^' [1 - vl V] (10) 

Since l+e^'^^^^^^g*^^^! vanishes due to the definition of the Dirac cones, the coupling 
term simplifies to 

C± = e-'^^'^'^'vW - e^^±^"^^^lV = ^[dy - ±iS^]. (11) 

Introducing a momentum operator q = — iV, our coupling term becomes C^ ~ 
fJ^Qx + '^Qy and the Hamiltonian spatial density is therefore: 

h.{r) = -hX{r) (12) 
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which is the expression of the massless Dirac Hamiltonian with Fermi velocity 
c = \/3Jd/2. Note that changing flavor is equivalent to changing the sign of Qx 
in the Hamiltonian. 

3. Extending the Dirac Hamiltonian by modifying the lattice parameters 

We now generalize the free Dirac Hamiltonian to include a variety of external fields [20] 

Hd = cd-p^ (3mc^ + (3Vcov (13) 

Here m is the mass of the particles and Vcov is the most general covariant potential 
containing scalar, vector, matrix, pseudoscalar, pseudovector and pseudotensor fields. 
The Dirac matrices a = 7^7 and ^5 = 7^, are defined in terms of the generators of 
the Clifford group, 7^. In 2+1 dimensions this is a set of 2 x 2 matrices satisfying 
the anticommutation relation {7^,7^} = 2r]^^ with r]^^ = (im^(l, — 1, — 1), where 
/i, z/ = 0, 1, 2. An adequate choice is a = a = {ax^ (Jy) and /3 = a^, which corresponds 
precisely to the Hamiltonian in equation ([g]). The Dirac Hamiltonian has the form 

Hd = ca-(p - eA{r)\^ eA^{r)l^ {mc^^V{f))(j^, (14) 

where A^{f) and A{r) are the usual scalar and vector potential that give rise to 
observable electric and magnetic fields, and V{r) is a scalar potential that mimics the 
effect of an imposed mass; e is an effective charge. Let us now show how to recover 



all terms in equation ( 14 ) by slightly perturbing the tight-binding model. 



3.1. Generating a mass term 

The simplest term that we can add to Hamiltonian equation ^ is an energy difference 
between atoms in sublattices A and B: 

^H = Y^ (I^L^m - |C^m) , (15) 

m 

which in the continuum limit becomes 

^^ = E / {^^U^)^a,{rl - |^J^(r)^5M(r)) d^r (16) 

or equivalent ly 

Sh{r)^a/-. (17) 



Comparing this with equation (14), we see that the position-independent term 



proportional to a^ can be identified with an effective mass. 

3.2. Abelian potentials 

We can extend our model by starting from equation Q and smoothly modifying the 
hoppings as J^ = J + Ci^rn where |e| <C | J|: 

SH = ^ ei^rn(^lj)m+vi + H.C. (18) 

m 

In the continuum limit this renders 

^^ = E/e^'^"'^'(^"*)*a^(^1*''^(^^+^^)d'r- (19) 
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In the notation of the previous section, this is equivalent to a change in the "couphng 
term" between pseudospins: 

therefore ahowing us to formahy derive an Abehan external potential 

H = cB-(q-A{r)\ (21) 

where A= {iiRe[5C^],lm[5C^]). 

3.3. Scalar fields 

Consider now perturbations which are diagonal in the internal space to the tight- 
binding Hamiltonian, such as those given by intralattice hoppings in equation (IT]): 

^H = ^ TAol^am+vi + TBhlj)m+vi + h.c. (22) 



Following the continuum limit performed in (3.1 ) the Hamiltonian in momentum space 
around the K point can then be written as 

^..,,.,-_(E<.-"T.(.l J.\. ,23) 

or purposefully rewritten (reabsorbing the constant ^^ ^-^KVi factor) as 

H' = ea ■ ,-- (lALl±Ij!m, _ (im_Ij!m,^ (24) 

which has both a term proportional to I representing an electric potential, (j){r) = 
^[r^(r) + r^(r)], and a term proportional to az which contributes to the effective 
mass. 

3.4' Flavor- coupling perturbations 

We have seen that spatial variations of the hopping elements emerge as an Abelian 
external field. These variations have to be small^ but it should be stated that, in order 
to keep flavors decoupled, they also have to be slowly varying. Otherwise, if there is a 
modulation with wavevector comparable to the order (i?+ — i?_), the flavor- coupling 
terms jj. ^ u in equation (^ do not vanish. What we suggest now is to introduce small 
perturbations whose wavelength is comparable to the lattice constant and thus bridge 
the difference in momentum between cones. For instance 

SH = Y1 ^Xx,i,m cos(r(^+ - K_))albm+y^ (25) 

m,i 

with a constant coupling strength Xx,i,m = Xx- These rapid oscillations only allow 
the survival of terms that couple different cones, cancelling all terms inside the same 
cone. The most general perturbation of this sort is 

6h{r) = xx(r-) (*l+*6- + *L*6+ + H.c.) , (26) 

where ^^ Xcc,*^"*^^^ = Xx is the spatial dependence of the slow envelope that 



surrounds our perturbation (25). Actually, this envelope can be "remodulated" in 



order to make flavor coupling also weakly spatial-dependent. Moreover, we can also 



introduce in equation (25) sine terms, e(r) = x?/ sin (r(i^+ — i^-)), which make the 



flavor-coupling term complex. 
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3.5. The complete Hamiltonian 

Since there are no pseudo-potentials (7^ = 1 in 2+1 dimensions), our idea of perturbing 
the tight-binding model parameters allows us to reconstruct all possible external 
potentials of the Dirac equation, plus an additional coupling between different types 
of particles. With all these elements we have the following effective single-particle 
Hamiltonian: 



ca{q- A{f)) + mazC^ + ^{f) XxC^x + Xy^y 

-c{a{q- A{r))y + mc'^a^ + (/)(r) 



h(f) = \ ^^w "^v^;; ^ ^^^^^^ ^y^yj xx^x^xy^y /27^ 



4. Trapping of atoms in an optical honeycomb lattice 

In this section we introduce optical lattices and explain how to construct a honeycomb 
lattice with two state-dependent triangular optical lattices. We discuss how to 
implement a tight-binding Hamiltonian and propose an experimental setup to obtain 
the Dirac field and the external fields using solely Raman and detuned lasers. 

4.1. From two triangular sub-lattices to a honeycomb lattice 

Recent advances in the development of high-aperture objectives and their integration 
in optical traps open the door to the generation of almost arbitrary two-dimensional 
potential landscapes for ultracold atoms. The basic idea is that off-resonant light 
may be used to tightly confine atoms in the maxima or minima of intensity, 
recreating sophisticated lattice models [21 . While until now those minima and 
maxima were generated through the interference of multiple laser beams [71 [21], a 
novel paradigm consists on shaping and organizing those intensity profiles by simply 
projecting sophisticated images on the two-dimensional focal plane of a lens. The 
first experiments along this line have reproduced the usual square lattice quantum 
simulations [22] and also demonstrated the first triangular lattices [23^ . 

In this work we are particularly interested in the last of those setups, which 
combines two triangular lattices [23] in the same plane. In this experiment the 
trapping laser beams are first diffracted by a holographic mask with a triangular 
pattern, selecting the first diffraction orders, which are then collected by a powerful 
lens to create the imprinted intensity pattern at its focal plane. The relative phases of 
the three diffracted beams may be independently controlled in a way that allows the 
displacement of the resulting triangular lattice (cf. figure [3|. When this procedure is 
applied to two laser beams that differ in frequency or polarization, it becomes possible 
to produce two triangular lattices, A and 5, which coexist on the same plane and have 
a tunable relative separation. The result is the original setup introduced in figure [l] 
in an abstract way, where now A and B are physically implemented by an optical 
potential. 

In order to jump from two independent triangular lattices to the setup introduced 
in figure [l] we need means to introduce the couplings between different lattice sites, 
that is hoppings from one lattice to another, or within the same lattice. The most 
direct way to implement this in an optical system is to use the two lattices A and B to 
trap atoms of the same species but in different internal states, \a) and \b) respectively. 
We also need a way to rotate between \a) and |6), which can be implemented via 
Raman transitions. In a situation with all these ingredients two kinds of hoppings 
are allowed for the atoms: on the one hand an atom in sub-lattice A {B) can tunnel 
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Figure 3. Proposal of experimental setup. DE is the diffraction element, LI and 
L2 are the required lenses to make the beams parallel and to focalize them on the 
lattice plane, Rl and R2 are the lasers used for the Raman transitions. Pi and 
P2 are phase plates used to change the relative phases, e^'^i'^, of the (T+ and a~ 
beams. 



between positions in its own sub-lattice, as contemplated in our model by the hopping 
parameter Ta (^b)- This parameter can be controlled by increasing or decreasing the 
intensity I a (Ib) of the laser beam generating each sublattice. On the other hand 
a Raman-induced change in the internal state of an atom from \a) to \b) (from \b) 
to I a)) will make the atom shift from sublattice A to sublattice B {B to A). These 
nearest-neighbour jumps can be different in each of the three possible directions, and 
are modelled by the hopping parameters Jo, Ji and J2 which are associated with the 
vectors vq , Vi and V2 shown in figure [l] 



4.2. Realization of two state- dependent triangular sub-lattices 

How do we implement in practice the ideas from the previous subsection, and in 
particular the coupling between lattices? The engineering and control of state- 
dependent lattices is a mature technology [24l [21] , which nevertheless requires some 
careful control of the atomic states and decoherence. We will briefiy describe how this 
works for fermionic alkaline atoms and how this integrates with the projected lattices 
scheme. 

As sketched in figure |4J it is possible to find a wavelength Xi falling between 
the Dl and D2 lines for which polarised light only traps atoms in one of the ground 
states manifolds. More precisely, the fine- structure energy levels of alkaline atoms 
are denoted by |L, J, ttij), where L is the electron angular momentum quantum 
number, J the combined orbital and spin momentum and rrij is its projection along 
the quantization axis. When we illuminate with cr+ polarised light at a frequency 
uj = ^{udi + ^D2), the ac-Stark shift that it induces on the |0, ^, — ^) cancels due 
to the positive and negative contribution of off-resonant Dl and D2 transitions. 
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Figure 4. Fine-structure energy levels of an alkaline atom. Not to scale. \di s-nd 
\d2 are the wavelengths of the Dl and D2 lines respectively. A^ is the wavelength 
of the optical lattice lasers. Dashed lines mark the cancelled ac-Stark effects: a 
cr+ (cr~ ) polarised photon incident on a |L = 0,J= |, mj = —\) (|05^5+^)) 



state would see both the 1 1, ^, ^) and the 1 1, |, I) (|l, |,-|)and |l, §,-§)) 
states in such a way that the ground level energy displacements cancel with each 
other. The net effect is then only due to (7~ (cr+), marked with continuous lines. 



C2 \J = \, 

we select 



The result is that circularly polarized (j^{(j ) light can only trap atoms in the 

|0'i+^> (|0'i-i)) states. 

In practice, however, the situation is more subtle because atoms also have some 
hyperfine structure, induced by the coupling between the electronic and nuclear 
angular momenta. Let us focus on the fermionic species ^Li, in line with previous 
proposals [23]. Out of the hyperfine ground-states |F, m^?) = ci | J = \^mj = ^) + 
— ^) where F = I-\-J being / the nuclear angular momentum of ^Li (/ = 1), 

^\ - |i -i\ - -.f^\^ l\^J-U -i\ ^nr\ l/)\-l^_3\_i|l_l\ 
"'/~|2' 2/ ~ Y3l2'2/^y3l2' 2 / ^^^^ \^/ ~ \ 2^ 2 / ~ "^ I 2 ' 2/' 

If V-^ and a V- are the intensity distributions that result by illuminating with light 
in the a~^ and a~ polarisations, the states \a) and \b) will feel the ac-Stark potentials 
Va = I V+ + I VI and Vb = V-. As shown in figure pi the combined potentials can 
look like two displaced triangular lattices by choosing the adequate relative phases 
between the diffracted beams. Furthermore, the relative depths can be controlled by 
changing the ratio V+/V_, but this may require a new tuning of the relative phases. 
The last thing is to induce state rotations between \a) and |6), which can be 
done via Raman transitions. So summarizing we have two triangular state-dependent 
sub-lattices A and B with intralattice hopping parameters Ta and F^. These two 
sub-lattices, combined with internal state rotations, give rise to a hexagonal lattice 
with interlattice hopping parameters J^. A scheme of a possible setup with all these 
elements is shown in figure |3] 



4.3. Perturbing the hopping parameters 

The previous discussion shows a feasible setup to obtain a honeycomb tight-binding 
Hamiltonian with neutral atoms in an optical lattice. In this section we take the 
scheme one step further so that the hopping parameters can be locally changed to 
incorporate mass terms and pseudofields to the Dirac Hamiltonian. We propose to 
optically control the Dirac field physics much in the way mechanical strains have been 
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Figure 5. Left: sub-lattice generated by Va 
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Right: sub-lattice 



generated by Vb = V_. Top: trapping potentials as functions of x and y, in 
arbitrary units. The red zones represent the potential minima where the atoms 
are trapped. It is apparent that both sub-lattices are triangular and displaced 
one with respect to the other. By superimposing both images one obtains the 
structure of an hexagonal lattice. Bottom: transversal cut corresponding to the 
white dashed line along the x direction. The difference in potential depth can be 
controlled by changing the ratio V^/V-^ but this will affect the relative positions 
of the potential minima. 



suggested in graphene sheets ^15l [161 [II]- Two proposals to implement variations 
are presented, the first of which is best in line with the setup shown in the previous 
section, the other being a more formal approach to general hopping calculations across 
potential wells. 

4.3.1. Change in Raman intensity Interlattice hoppings are assisted by an external 
laser beam inducing Raman transitions between the internal states of the atoms, 
thus creating a transition amplitude which scales with the intensity of the beam. 
By changing the spatial dependence of this amplitude one can overprint a first- 
neighbours hopping perturbation ei{r). This method is the easiest to implement in this 
setup, but is very limited and unable to act on intralattice (next-to-nearest-neighbour) 
transitions. Therefore it is the chosen method to add Abelian vector potentials and 
non-Abelian coupling terms. 



4.3.2. Lattice distortions A more general analysis can be made by considering spatial 
distortions of the lattice, that is, either through a relative displacement of one site 
with respect to its neighbor or by adjusting width of the individual potential wells. 
The displacement may by achieved by controlling the relative phases of the conforming 
beams through the phase modulators PI and P2 in figurejs] The width of the potential 
wells can be modified by a change in the height of the confining potential, using a mask 
which modulates the intensity of the laser beams. We now develop a simple model for 
calculating the influence of these two modifications in the hopping parameters. 

We assume that the ideal (unperturbed) hopping parameters between wells 
behave as if lattice sites were spatially separated harmonic oscillators in an original 
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reference frame r ^ The lattice distortion described locally by the transformation 
f ' = Ar renders the following single- well Hamiltonian: 

H= A+ r^A^Ar (28) 

2m 2 ^ ^ 

where the frequency of the site trap is also changed uj -^ cD(r), due to modulations 
of the intensity of the confining beam (the kinetic energy terms are not transformed 
since they are expressed in "real" space, while the potential term is meant to look like 
a perturbed harmonic oscillator). The localized Wannier function in the well is 

^Gs{x) = -^ exp(-^r^5r) (29) 

where B = VA^, a^ = -t and TV = ^^. Notice that B is positive-definite by 

^^'^ \/\B\ 

construction and a is r-dependent. If we assume no vibrational levels can be excited, 
the hopping parameter between first-neighbour wells separated by a vector a is simply 
proportional to the overlap between wavefunctions 

J ~ Jd^rrosi^^Gsir- a) = exp(-^). (30) 

This will be valid as long as variations are smooth and small, corresponding to 
the wavefunctions being good approximations to the ground state of the potential 
well and the restriction to the lowest vibrational level. Therefore the dislocation 
must be smah: Ac^I^eC^G = A^A c^ I + e(C + C^), defining C as our 
generator. Let us call the unperturbed vector connecting both wells vq^ so that the 
actual vector becomes a = A~^vo = (I — eC)vo. The bilinear matrix expands as 
B = VG c^ ^/ITe(cT~C^ ::^ I + ^e(C + C^). Substituting these values yields: 

J ^ exp( :r^r) ^ J = Joexp(-eO), 

6=-^v^{C + C^))vo. (31) 

This result can be readily interpreted. There are two contributions: one {C + C^) 
which modifies the hopping due to the change in distance between wells, and another 
{oj{r)) which refers to the intensity profile of the confining laser beam. We can therefore 
tune these parameters, e.g. displacing or rotating one sublattice on top of the other, 
to allow for variations of the tight-binding hoppings. As explained before, all these 
changes are used to simulate Abelian and non-Abelian external gauge fields. 

5. Experimental detection 

The most direct consequence of the appearance of Abelian and non-Abelian external 
fields is the distortion of the energy bands, ranging from the movement of the Dirac 
cones in momentum space to the appearance of a gap in the spectrum. Two methods 
are proposed here to observe these changes: the first one is the measurement of the 
momenta after removing the confining potential so as to explore the form of the energy 
bands; the second one focuses on measuring how these fields affect the dynamics of a 
single particle (or small group of particles) moving in the lattice. 
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Figure 6. Expected momenta distribution for B (left) and A (right) atoms 
associated with the E- and E-^ energy bands respectively (see text). The two 
top images correspond to the unperturbed honeycomb tight-binding Hamiltonian, 
and the bottom ones correspond to the Abelian potential case described in section 
[3] due to a perturbation in Jq. 

5.1. Time of flight images 

It is possible to probe the atomic population within the first Brillouin zone (BZ) of 
an optical lattice. The procedure consists of the following steps: (i) adiabatically 
switching off the lattice potentials so that the atoms quasimomenta are converted into 
real momenta, (ii) then letting the atoms expand freely during a certain time of fiight 
and (iii) finally taking an absorption image of the expanded cloud [25 . 

In our setup the situation is a bit more complicated. Let us denote by |+,/c) 
and | — ,/c) the eigenstates associated with upper and lower energy bands, E_^{k) and 
E_{k), of the tight-binding Hamiltonian (l2|. In a ground state with a Fermi energy 
slightly above zero, we would have many atoms in the lower band sector occupying the 
whole of the first BZ and only few atoms in the upper band, concentrated around Dirac 
points (See figure [6|i). Despite the difference, by means of absorption images one can 
not distinguish between |+) and |— ) states as both have components corresponding 
to the |a),|6) internal states of the atoms. 

In order to picture the Dirac cones we need a method that discriminates between 
energy bands, say, transforming all the |+) states into |a)'s and all the |— ) states into 
|6)'s, while preserving the momentum, k. The adiabatic theorem provides us the way 
for doing this. Our starting point is the honeycomb lattice Hamiltonian in momentum 
space 

^ = -^ E 4 [/(^)'^^ + /*(^)'^"] H' (32) 

k 

written with the pseudospin structure u^ = (a^, 6^), the couplings f{k) = 1 + e*^'^^ + 
^-ik-v2 ^ ^^^ ^Y^Q Pauli ladder operators a^. Note how the Hamiltonian is a composition 
of commuting terms for each value of the momentum, k. We will adiabatically distort 
all terms, following the route in figure [7| which consists on first adding a mass term. 
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Figure 7. Adiabatic path for a transformation of the |+}, |— } states into \a), 
\b) states. The hamiltonian of the system is written as if ~ s(J,m) • a, so 
^initial = s(J,0) lies in the x,y plane while s final = s(0,m) lies along the z 
axis. During the first step m is adiabatically increased, and during the second J 
is adiabatically decreased until zero. 



Table 1. summary table of all the considered perturbations of the hopping 
parameters with their physical effects. 



Lattice modification 


Hamiltonian modification 


Effect 


Ji^J^a i=0,l,2 
small hopping parameter 
modification 


A = A(cos(K • Vi) , sin(i? • ^T^)) 


Abelian field 


energy difference between 
A and B 


H ^ a ■ q-\- ^CFz 


Mass term 


intrasublattice triangular 
hoppings 


H^a-q-%{To. + Tp)l -%{T^- Tp)a, 


Scalar field 


e(^) = Xx cos(r(i^+ - ^-)) 


6h{f) = x4^i+^b- + ^I_^6+ + H.c.) 


Flavor 
coupling 



~ 771(7^, and then decreasing J down to zero. The protocol maps the two eigenstates 
of the initial Hamiltonian, |+) and |— ), to \a) and |6), accurately. 

This technique allow us to experimentally observe the effects of an Abelian 
potential described in section [3) as represented in figure [gJd, where the Dirac cone 
displacement manifests as a deformation of the |a)-momentum distribution. In general 
all other effects, which are summarized in Table [l] could be observed by similar 
methods. 



5.2. Few particle dynamics 

The method in the previous subsection gives us access, among other things, to the 
best known experimental observable in condensed matter physics, which is the density 
of states. It has been the easiest to measure and therefore has become the default 
choice in optical lattice simulations of solid state physics. However, the distinctive 
characteristics of the excitations in a half-filled hexagonal lattice may be worth an 
extra effort: trying to observe the behaviour of a group of particles with well-defined 
momentum obeying a complete Dirac equation will reveal some of the distinctive 
features of the simulated fields. 



Simulating the Dirac field with cold atoms in a honeycomb lattice 



15 



a) 



Pe 




Figure 8. a) Excitation probability for three different cases: No gap (m=0, 
solid line), medium gap (m=V, dashed line) and large gap (m=2.5V, dotted line) 
against time. While the gapless state gets immediately promoted once it reaches 
the Dirac points, the gapped states have reduced positive energy contributions, 
b) Change in {(7x)'. {(Jx{t = t)) — {(Jx{t = —r)) depending on the value of the 
effective mass. We compare the change in {ax) by an exact simulation of the 
particle dynamics with the Landau-Zener formula (purple). 



As an example we suggest using the Klein tunneling effect [26 to probe and 
measure the energy gaps between the two bands. Let us add a uniform electric field 
pointing along one direction, say x, described by a linearly growing potential, Vx. As 
explained in Ref. |27j, the effect of this potential is to accelerate particles, continuously 
increasing their momenta in time, kx{t) = kx{0) — Vt/h^ until the particle reaches a 
boundary of the Brillouin zone. At this point two things may happen. If the particle 
is far away from a Dirac singularity, or the minimum gap between energy bands (the 
effective mass m) is large compared to the acceleration, mc'^ ^ V, the particle will 
simply reappear through the opposite side of the Brillouin zone, reversing its velocity 
and performing the so called Bloch oscillations. However, if the particle hits against 
the proximities of the K± points and the gap is small, mc'^ <C V/h^ the particle will 
experience a Landau-Zener process and jump to the opposite energy band, maintaining 
its group velocity. 

As shown in figure [8J we have simulated numerically this process for different 
masses of the Dirac field. As a signature of the jump between bands we simply use 
the expectation value of k-a, which is directly correlated to band excitation probability 
but easier to measure. Note how for the gapless phase, m = 0, there is a perfect jump, 
and how the probability is well approximated by the Landau-Zener formula, which 
allows us to reverse-engineer the experiment and fit the value of m. 

How would this be implemented in an experiment? The idea would be to do the 
same optical lattice setup with a small number of fermions or bosons cooled down to 
the lowest value of the momenta, A: = 0, in a honeycomb lattice that implements the 
desired Dirac Hamiltonian, with or without effective mass. One would then activate 
an electric field along the direction w, for a certain time t. After this time one would 
measure the state of the atoms, or more precisely the expected value (a). By changing 
the duration of the field and its intensity, and monitoring the changes in /c • a, we will 
be able to reconstruct not only the Klein effect but also the whole spin texture of the 
bands. 
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6. Conclusions 

We have presented a system composed by two triangular sub-lattices which conform a 
honeycomb lattice. We have seen how the different terms of the two-dimensional Dirac 
Hamiltonian arise by perturbing the lattice hopping parameters. We have also made 
an experimental proposal to implement the model, and discussed how the effects of 
the different terms in the Dirac Hamiltonian can be observed by detecting either the 
momenta distribution of the atoms or the relative densities of atoms at each sublattice. 
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